import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib as mpl
import scienceplots
# Science plot package + Dracula theme
plt.style.use(['science', 'dark_background'])
plt.rcParams['axes.facecolor'] = '#282a36'
plt.rcParams['figure.facecolor'] = '#282a36'
colorcycle = ['#8be9fd', '#ff79c6', '#50fa7b', '#bd93f9', '#ffb86c', '#ff5555', '#f1fa8c',
'#6272a4']
plt.rcParams['axes.prop_cycle'] = mpl.cycler(color=colorcycle)
white = '#f8f8f2' # foreground
# change dpi
plt.rcParams['figure.dpi'] = 512
The 'method':
import cv2 as cv
# load image
img_1rgb = cv.imread('../images/hw3_images/emb1_rgb.tif')
img_1bw = cv.imread('../images/hw3_images/emb1_2c.tif')
# image dimensions
print(img_1rgb.shape) # (height, width, channels)
print(img_1bw.shape) # (height, width, channels
# convert top right corner 500x500 to black
img_1rgb[:500, -500:] = 0
img_1bw[:500, -500:] = 0
# Opencv uses BGR, convert to RGB
rgb_image = cv.cvtColor(img_1rgb, cv.COLOR_BGR2RGB)
# 2 channel image
# Define a threshold for identifying top and bottom bounds
threshold = 60 # Adjust as needed
# finding the top and bottom bounds for cell
top_bound = None
bottom_bound = None
for row in range(img_1bw.shape[0]):
for col in range(img_1bw.shape[1]):
if img_1bw[row, col, 0] > threshold:
if top_bound is None:
top_bound = row
bottom_bound = row
print(top_bound, bottom_bound)
# Crop the image to the bounds
cropped_image = rgb_image[top_bound:bottom_bound, :]
plt.imshow(cropped_image)
(2810, 2719, 3) (2810, 2719, 3) 163 2617
<matplotlib.image.AxesImage at 0x29a0d5d90>
print(img_1bw.shape)
print(img_1bw[334, 433, 0])
(2810, 2719, 3) 29
# # Array to store left and right bounds for each row
bounds = []
# Iterate through each row to find the left and right bounds
for row in range(top_bound, bottom_bound):
# Find left and right bounds for the current row
left_bound = None
right_bound = None
for col in range(img_1bw.shape[1]):
# check if white in left and right bounds
if img_1bw[row, col][0] > threshold:
if left_bound is None:
left_bound = col
for col in range(img_1bw.shape[1]-1, -1, -1):
if img_1bw[row, col][0] > threshold:
if right_bound is None:
right_bound = col
bounds.append((left_bound, right_bound))
# bounds to numpy array
bounds = np.array(bounds)
# print(bounds)
# delete first element
bounds = bounds[1:]
# print(bounds)
green_channel = rgb_image[:, :, 1]
# Calculate the average green value within the bounds for each row
average_green_values = []
i = 0
for row in range(top_bound, bottom_bound-1):
left_bound, right_bound = bounds[i]
average_green_values.append(np.mean(green_channel[row, left_bound:right_bound]))
i += 1
# Plot the average green values
plt.plot(average_green_values)
plt.show()
img_2bw = cv.imread('/Users/junseo/washu/spring2024/phys481/images/hw3_images/emb2_2c.tif')
img_2rgb = cv.imread('/Users/junseo/washu/spring2024/phys481/images/hw3_images/emb2_rbg.tif')
# image dimensions
print(img_2rgb.shape) # (height, width, channels)
print(img_2bw.shape) # (height, width, channels
# threshold for identifying top and bottom bounds
threshold2 = 100 # Adjust as needed
# finding the top and bottom bounds for cell
top_bound2 = None
bottom_bound2 = None
for row in range(img_2bw.shape[0]):
for col in range(img_2bw.shape[1]):
if img_2bw[row, col, 2] > threshold2:
if top_bound2 is None:
top_bound2 = row
bottom_bound2 = row
print(top_bound2, bottom_bound2)
# Crop the image to the bounds
cropped_image2 = img_2rgb[top_bound2:bottom_bound2, :]
plt.imshow(cropped_image2)
(3682, 1992, 3) (3682, 1992, 3) 197 3483
<matplotlib.image.AxesImage at 0x173b3bf50>
# draw the image for the 2 channels in grayscale
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(img_2bw[:, :, 1], cmap='gray')
ax[1].imshow(img_2bw[:, :, 2], cmap='gray')
plt.show()
# Array to store left and right bounds for each row
bounds2 = []
# Iterate through each row to find the left and right bounds
for row in range(top_bound2, bottom_bound2):
# Find left and right bounds for the current row
left_bound = None
right_bound = None
for col in range(img_2bw.shape[1]):
# check if white in left and right bounds
if img_2bw[row, col][2] > threshold2:
if left_bound is None:
left_bound = col
for col in range(img_2bw.shape[1]-1, -1, -1):
if img_2bw[row, col][2] > threshold2:
if right_bound is None:
right_bound = col
bounds2.append((left_bound, right_bound))
# bounds to numpy array
bounds2 = np.array(bounds2)
# print(bounds2)
# delete rows with no bounds
bounds2 = bounds2[np.where(bounds2[:, 0] != bounds2[:, 1])]
img_2rgb = cv.cvtColor(img_2rgb, cv.COLOR_BGR2RGB)
green_channel2 = img_2rgb[:, :, 1]
# Calculate the average green value within the bounds for each row
average_green_values2 = []
i = 0
for row in range(top_bound2, bottom_bound2-2):
left_bound2, right_bound2 = bounds2[i]
average_green_values2.append(np.mean(green_channel2[row, left_bound2:right_bound2]))
i += 1
# Plot the average green values
plt.plot(average_green_values2)
plt.show()
# error bars for emb1
std_dev1 = []
i_1 = 0
for row in range(top_bound, bottom_bound-1):
left_bound, right_bound = bounds[i_1]
std_dev1.append(np.std(green_channel[row, left_bound:right_bound]))
i_1 += 1
plt.errorbar(range(len(average_green_values)), average_green_values, yerr=std_dev1, ecolor=colorcycle[1], label='error bars')
plt.title('EMB 1: Average Green Value with Error Bars')
plt.legend()
plt.show()
# print(bounds2)
# calculating error bars for emb 2
# standard deviation
std_dev2 = []
i_2 = 0
for row in range(top_bound2, bottom_bound2-2):
left_bound2, right_bound2 = bounds2[i_2]
std_dev2.append(np.std(green_channel2[row, left_bound2:right_bound2]))
i_2 += 1
# print(std_dev2)
# Plot the average green values with error bars at x = 1000, 1500, 2000, 2500, 3000
plt.errorbar(range(len(average_green_values2)), average_green_values2, yerr=std_dev2, ecolor=colorcycle[1], label='error bars')
plt.title('EMB 2: Average Green Value with Error Bars')
plt.legend()
plt.show()
The issues/ shortcomings:
This issue is amplified in embryo 1, where the embryo lies on an oblique angle, and we can see this within the error bars in the interval close to the middle of the embryo. We could improve upon this by either rotating the image so that the embryo is parallel to the $y$ axis, or by using a more sophisticated method where we transform the values on to the embryo's long body axis.
inaccuracy of the bounds for the cells body.
for each embryo. This inaccurately finds the bounds of the embryo which may have also resulted in the error bars being larger than they should be. We could improve this by applying a better filter for noise and an edge detection algorithm (e.g. Canny edge detection) to find the bounds of the embryo.
We need to somehow match the $x$ values of the two embryos, and since embryo 2 is a larger array, we will have to guess/interpolate what the $x$ values of embryo 1 would be if it were the same size.
After we have 'normalized' the $x$ axis of the two embryos, we could possible compare the models by taking the difference of each $y$ value at each row along the long body axis. We would expect that this difference would be close to zero if the models are similar, and the model kind of shows this relation (although the error bars are quite large).
# turn into numpy arrays
average_green_values = np.array(average_green_values)
average_green_values2 = np.array(average_green_values2)
# turn into numpy arrays
average_green_values = np.array(average_green_values)
average_green_values2 = np.array(average_green_values2)
# determine the length of the longer array
max_len = len(average_green_values2)
# create a new x array with the length of the longer array
x = np.linspace(0, 1, max_len)
# create new x arrays for the original data
x1 = np.linspace(0, 1, len(average_green_values))
x2 = np.linspace(0, 1, len(average_green_values2))
# interpolate the data to the new x array
ynew = np.interp(x, x1, average_green_values)
ynew2 = np.interp(x, x2, average_green_values2)
# plot the interpolated data
plt.plot(x, ynew, label='Embryo 1')
plt.plot(x, ynew2, label='Embryo 2')
plt.title('Interpolated Data for EMB 1 and EMB 2')
plt.legend()
plt.show()
# Comparing profiles of EMB 1 and EMB 2 by taking the difference of mean values
diff = ynew2 - ynew
# Standard deviation of the difference
std_dev_diff = np.std(diff)
error_x = [500, 1000, 1500, 2000, 2500, 3000]
error_y = [diff[500], diff[1000], diff[1500], diff[2000], diff[2500], diff[3000]]
plt.plot(diff)
plt.title('Difference of Mean Values of EMB 1 and EMB 2')
plt.errorbar(error_x, error_y, yerr=std_dev_diff, ecolor=colorcycle[1], fmt='none')
plt.legend()
plt.show()
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.